double energiaLinkCellList(caja *bx)
{
	//Si calcStress = true, se calcula la parte correspondiente al potencial
	//del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
	int c, c1;
	int iii,jjj;
	int shift[3];
	int count=0;//Numero de pares interactuantes
	double dist,dist2,dist3,potdm, ePoten;
	double Co,Co2;
	double rc = bx->rcut;
	double rc2=rc*rc;
	double virial=0;
	double f2;

	for(int i=0;i<bx->npart;i++)
	{
		bx->part[i].energia=0.0;
	}

	//Inicializamos la energia potencial
	ePoten=0.;
	
	//Recorremos todas las celdas
	for (int ix=0;ix<bx->ncell;ix++)
	{
		for (int iy=0;iy<bx->ncell;iy++)
		{
			for (int iz=0;iz<bx->ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*bx->ncell2+iy*bx->ncell+iz);

				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+bx->ncell);
					}
					else
					{
						if (ixred>=bx->ncell)
						{
							shift[0] = int(ixred-bx->ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+bx->ncell);
						}
						else
						{
							if(iyred>=bx->ncell)
							{
								shift[1] = int(iyred-bx->ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+bx->ncell);
							}
							else
							{
								if(izred>=bx->ncell)
								{
									shift[2] = int(izred-bx->ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell
								+((shift[2]+bx->ncell)%bx->ncell));
							///Buscamos la primera particula iii de la celda c
							iii = bx->head[c];

							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =bx->head[c1];
								
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
					
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										///Calculamos la distancia con condiciones periodicas entre iii y jjj
										dist2= bx->distancia2(iii,jjj);
										//Si la distancia es menor que el cutoff calculamos la interaccion
										if (dist2 < rc2)
										{
											dist=sqrt(dist2);
											dist3=dist*dist*dist;
											count++;
											if (dist<4.0)
											{
												//Para el potencial
												f2=0.00375*(bx->u0*exp(-50.*(dist-2.)));
											}
											else
											{
												f2=0.0;
											}

											///Calculamos la fuerza entre iii y jjj
											Co=bx->dr[2]/dist;
											Co2=Co*Co;				

											potdm = f2+(bx->u0)*(1.0-3*Co2)/(dist3);

		

											bx->part[iii].energia+=potdm;
											bx->part[jjj].energia+=potdm;
											//Acualizamos la energia potencial
											ePoten+=(potdm);
										}
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = bx->link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = bx->link[iii];
							}
						}
					}
				}
			}
		}
	}

	//Introducimos la correccion del cutoff del potencial
	bx->epoten=ePoten/bx->npart;
	return ePoten;
}

double calcEnergiaParticle(caja *bx, int index)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	double dist, dist2, dist3, potdm, ePoten;
	double rcut2=bx->rcut*bx->rcut;
	double f2, Co, Co2;
	int count=0;

	ePoten=0.;

	c = bx->part[index].celda;
	///Recorremos las particulas de esta celda y las celdas vecinas
	for (int ixred=bx->part[index].celdaxyz[0]-1;ixred<=bx->part[index].celdaxyz[0]+1;ixred++)
	{
		/*shift: condiciones de contorno periodicas tenidas en cuenta para
		localizar las celdas vecinas*/
		if (ixred < 0)
		{
			shift[0] = int(ixred+bx->ncell);
		}
		else
		{
			if (ixred>=bx->ncell)
			{
				shift[0] = int(ixred-bx->ncell); 
			}
			else
			{
				shift[0] = int(ixred);
			}
		}
		for (int iyred=bx->part[index].celdaxyz[1]-1;iyred<=bx->part[index].celdaxyz[1]+1;iyred++)
		{
			if (iyred < 0)
			{
				shift[1] = int(iyred+bx->ncell);
			}
			else
			{
				if(iyred>=bx->ncell)
				{
					shift[1] = int(iyred-bx->ncell);
				}
				else
				{
					shift[1] = int(iyred);
				}
			}
			for (int izred=bx->part[index].celdaxyz[2]-1;izred<=bx->part[index].celdaxyz[2]+1;izred++)
			{
				if (izred < 0)
				{
					shift[2] = int(izred+bx->ncell);
				}
				else
				{
					if(izred>=bx->ncell)
					{
						shift[2] = int(izred-bx->ncell);
					}
					else
					{
						shift[2] = int(izred) ;
					}
				}
				///Calculamos el indice c1 de la celda vecina
				c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell+((shift[2]+bx->ncell)%bx->ncell));
				///El indice index ahora es iii
				iii = index;

				///Buscamos la primera particula jjj de la celda c1
				jjj =bx->head[c1];

				while (jjj != -1) ///jjj=-1->No hay particulas en la celda
				{
					if(iii<jjj)
					{

						///Calculamos la distancia con condiciones periodicas entre iii y jjj
						dist2= bx->distancia2(iii,jjj);
						//Si la distancia es menor que el cutoff calculamos la interaccion
						if (dist2 < rcut2 )
						{
							dist=sqrt(dist2);
							dist3=dist*dist*dist;
							count++;
							if (dist<4.0)
							{
								//Para el potencial
								f2=-0.00375*bx->u0*exp(-50.*(dist-2.));
							}
							else
							{
								f2=0.0;
							}

							///Calculamos la fuerza entre iii y jjj
							Co=bx->dr[2]/dist;
							Co2=Co*Co;																	

							potdm = f2+(bx->u0)*(1.0-3*Co2)/(dist3);
							//Acualizamos la energia potencial
							ePoten+=potdm;							
						}
					}
						//Pasamos a la siguiente particula de la celda c1
						jjj = bx->link[jjj];
				}
			}
		}
	}
	return ePoten;
}
